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We investigate the possibility that the surface of a strongly coupled electron-phonon system behaves differ- 
ently from the bulk when the relevant parameters are inhomogeneous due to the presence of the interface. We 
consider parameter variations which make the surface either more metallic or more insulating than the bulk. 
While it appears impossible to stabilize a truly insulating surface when the bulk is metallic, the opposite situ- 
ation can be realized. A metallic surface can indeed be decoupled from a bipolaronic insulator realized in the 
bulk. 

PACS numbers: 71.38.-k, 71.30.+h, 73.20.-r, 71.38.Ht 



I. INTRODUCTION 

The interest in the electronic properties of surfaces and in- 
terfaces is growing due to the increasing ability to engineer 
interfaces between correlated materials and to accurately mea- 
sure surface and bulk properties. A number of discrepancies 
have been reported between bulk and surface properties of 
complex materials 1 - 2 , while interfaces between different ma- 
terials can lead to surprising properties. A notable example is 
the metallic interface between the two insulators LaTiC>3 and 
SrTi0 3 3 . 

On the theoretical side, the investigation of the effects of 
surfaces and interfaces has been focused on Hubbard-type 
models, in which local repulsion correlates the electronic 
motion eventually leads, for commensurate densities, to a 
Mott insulating state when the Coulomb interaction is suffi- 
ciently large. These studies have either used extensions 4 ' 5 of 
the dynamical mean-field theory (DMFT) 6 , a theoretical ap- 
proach which has provided the first unified scenario of the 
Mott transition, or variational approaches. 7 Studies of solid- 
vacuum interfaces have unveiled the possibility of surface 
ferromagnetism 8 , and have described the penetration depth of 
a bulk metallic phase into an otherwise insulating surface. 9 In- 
deed, Borghi et al? have shown the existence of a dead layer, 
due to an exponential penetration of metallic excitations. 

Another localizing effect which affects the properties of 
electrons in solids is the electron-phonon (e-ph) interaction. 
Also in this case quantum fluctuations inherent to the low 
dimensionality of surfaces and interfaces and strong interac- 
tions can stabilize novel ground states that are distinct from 
the bulk. In Ref . 1 the freezing of a bulk phonon at the sur- 
face has been invoked as the source of remarkable electronic 
properties. 

Similarly to the case of repulsive electron-electron inter- 
actions, important insights into the problem of strongly cou- 
pled e-ph systems have been gained by DMFT 6 . Studies 
of the Holstein model in a homogeneous bulk system using 
DMFT 10-12 show that as the e-ph interaction increases, the 
conduction electrons progressively lose their mobility, even- 
tually evolving into a polaronic state in which the presence of 
an electron is associated with a finite lattice distortion. The 



same e-ph coupling can cause any two polarons to attract and 
form a bound pair in real space, called a bipolaron. When the 
number of carriers equals the number of sites, i.e., the lattice is 
half-filled, bipolaron formation causes the system to undergo 
a continuous (at zero temperature) metal to insulator transition 
at a critical e-ph coupling. 

We have investigated the effect of a solid-vacuum inter- 
face on this scenario in a previous paper, 13 in which, in order 
to focus on the purely geometrical aspect of the problem we 
have considered the same parameters in the surface and in the 
bulk. Even for uniform parameters, the band narrowing at the 
surface 14 causes a reduction in quasiparticle weight relative 
to the bulk, i.e., the surface is less metallic. Therefore, upon 
increasing the e-ph interaction strength the polaron crossover 
takes place first on the surface layer. Nonetheless, for uniform 
model parameters, enhanced correlation effects at the surface 
are not sufficient to turn the surface insulating before the bulk 
(i.e., for a smaller e-ph coupling) and a single metal-insulator 
transition occurs at the critical coupling for the infinite system 

13 

9c — 9c,bulk- 

Besides the geometrical effect of missing neighbors, the 
surface properties are complicated by the fact that the micro- 
scopic interactions close to the surface have a value which 
may differ significantly from that in the bulk. In the Holstein 
model a modification of the hopping as well as of the e-ph 
coupling strength in the vicinity of the surface should be ex- 
pected for any real systems. The relaxation of the interlayer 
distance, for example, can cause an enhancement or decrease 
in the hopping integrals at the surface. In this work we shall 
extend the analysis of Ref. 13 to nonuniform model parameters 
and investigate the possibility of the occurrence of a metallic 
surface concurrent with a bulk bipolaronic insulator or of a 
bipolaronic insulating surface concurrent with a normal metal 
in the bulk. We will consider the half-filled case, in which an 
actual phase transition can be observed, even if particle-hole 
symmetry forbids charge transfer between surface and bulk. 

The paper is organized as follows. In Sec. II we intro- 
duce the model Hamiltonian, which is a semi-infinite Hol- 
stein model with layer dependent parameters. In addition, we 
briefly describe the embedding approach for DMFT. Results 
for a range of modified surface parameters are presented and 
discussed in Sec. III. In Sec. IV we summarize with some 
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concluding remarks. 



II. MODEL AND METHOD 

We investigate the Holstein model on a three-dimensional, 
bipartite simple-cubic (sc) lattice with nearest-neighbor hop- 
ping. The lattice is cut along a plane perpendicular to one 
of the coordinate axes, e.g., the z-axis [sc(OOl) surface]. The 
system is considered to be built up by two-dimensional lay- 
ers parallel to the surface. Accordingly, the position vector 
to a particular site in the semi-infinite lattice is written as 
Rsite = fi + R a - Here R a stands for the coordinate ori- 
gin in the layer a and the layer index runs from a = 1 for the 
topmost surface layer to infinity, r; is the position vector with 
respect to a layer-dependent origin, and runs over the sites 
within the layer. Each lattice site is then labeled by indices i 
and a. In this notation, the Hamiltonian reads: 

{ia,jf3)a ia 



^g a {n la - 1) (b\ a + b ia j, 



(1) 



where c; QCT [c\ a ^j and 6 iQ f^Lj are ' respectively, destruc- 
tion (creation) operators for electrons with spin a and local 
vibrons of frequency Hq on site i of the a layer. The electron 
density on site ia is denoted n; Q , ti a ,^ is the hopping matrix 
element between two nearest-neighbor sites, and g a denotes 
the layer-dependent electron-phonon coupling strength. We 
fix the energy scale by setting tu a jp\ =t=l for a, /3 ^ 1. 

To solve our model, we use an extension of DMFT to 
inhomogeneous systems called the embedding approach for 
DMFT. 5 In this scheme, the layered structure is partitioned 
into a surface region which includes the first N layers, and 
the adjacent semi-infinite bulk region (substrate) which is cou- 
pled to it (see Fig. 1). The surface corresponds to the region 
where one expects different properties relative to the bulk. It 
is shown next that the influence of the semi-infinite substrate 
on the surface region can be described in terms of an energy- 
dependent embedding potential. This can be viewed as an 
additional self-energy due to the transitions between the sur- 
face and the substrate. Because of translational symmetry in 
the plane parallel to the interface, the embedding potential of 
the substrate is diagonal with respect to the two-dimensional 
wave vector k = (k x , k y ) and can be expressed as an N x N 
matrix. 

By defining A(k, ioj n ) = [(iu n + — e (k) — , 
the equation for the Green's function is given by: 



A(k, iw n )G(k, iu n ) 



While the surface region consists of only N layers, the matri- 
ces corresponding to the Green's function are infinite dimen- 
sional due to the semi-infinite substrate. In Eq. (2), H(iuj n ) is 
the self-energy matrix, which in the framework of single-site 
DMFT, is local [i.e., H(iu! n ) al 3 = S a ,(iw„)5 Q ,^] and indepen- 
dent of wave vectors, k. The chemical potential is given by 



SbuikM S bu i k (w) Sjv(w) Ei(w) 



tl2 



gi 



N+1 N 2 1 



substrate 



surface region 



FIG. 1: Geometry of the (001) surface of a simple cubic lattice. The 
index a (horizontal axis) labels the layers parallel to surface, a = 1 
refers to topmost layer. In the embedding approach for DMFT the 
system is divided into a surface region of N layers and a semi-infinite 
substrate. The physical properties in the substrate (e.g., electron self- 
energy) are not layer-dependent and coincide with those of the bulk. 



/i and e(k) is the two-dimensional dispersion relation, which 
includes information about the surface geometry. The e(k) 
matrix for a surface cutting a simple cubic lattice with a plane 
perpendicular to the z direction [sc(001) surface] assumes the 
following form: 4 



e(k) 



*ne||(k) ti a ei_(k) 

taie±(k) *aae||(k) t 23 e^(k) 

t 32 e ± (k) t33C||(k) ••• 





(3) 



The intralayer (parallel) hopping and the interlayer (perpen- 
dicular) hopping are specified by t aa eu(h) and t a pe±(k), 
respectively, 15 with 
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-2[cos(fc x ) +cos(k y )], |e_i_(k)| 



1. 



(4) 



Enforcing the separation between the surface (S) layers and 
the substrate (L), we can write Eq. (2) in a block form: 



All Als 
Asl A S s 



Gll Gls 
Gsl Gss 



1 
1 



(5) 



It should be noted that A LS 



Asl and Als is a sparse ma- 
trix independent of k and uj n . For nearest-neighbor hopping, 
Als nas only one nonzero entry corresponding to the hop- 
ping between the lowest layer of the surface and the top of the 
substrate. We can solve Eq. (5) to obtain the surface Green's 
function Gss- One obtains the following relation between 
N x N matrices defined on the surface block: 



ss - A S lA ll Als) G S s 



(6) 



The second term in the parenthesis of Eq. (6) defines the em- 
bedding potential due to coupling of the surface region to the 
(2) substrate. By definition A^ L is the Green's function of the 



substrate decoupled from the surface 

G(k,io; n )= [(iw n +/x)l - e(k) - S(w n )] _1 . (7) 
The embedding potential then reads 

S(k, iu) n ) = A SL G(k, iuj n )A LS - (8) 
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Since Als is nonzero only between nearest-neighbor layers 
of the substrate and surface regions, only the Green's func- 
tion of the first layer of the substrate, 14 i.e., the first entry 
Gn(k, iuj n ) of Eq. (7), is needed to calculate the embed- 
ding potential. Gn(k, i(J n ) is computable directly using a 
recursive relation. 14 The self-energy appearing in Eq. (7) is 
obtained through a standard DMFT calculation for the bulk 
crystal corresponding to the substrate. 

After constructing the embedding potential of the substrate, 
S(k, iu) n ), we can compute the self-energy of the surface 
layers by DMFT. This can be achieved via the following 
steps: (i) we associating an effective impurity model with each 
layer in the surface region, and solve them by using an im- 
purity solver to find the layer-dependent local self-energies, 
H a (iuj n ). Then we construct the surface region self-energy 
matrix which is diagonal in layer indices (a, (3) with the ele- 
ments, Y, a p{iu n ) = T, a (ioj n )S a j3, (ii) we calculate the on-site 
layer-dependent Green's function via the following relation: 

G a (iu> n ) = 

^ ((i^n +/")! - e(k) - S(k,ioj n ) - S(iw„)) QQ ' ^ 

where the N x N e(k) matrix is given by Eq. (3). (iii) We 
implement the DMFT self-consistency relation for each layer, 

Q%{ioj n ) = + Ti a (i(jj n )] , which determines the 

bath parameters for the new effective impurity model. These 
steps have to be repeated until self-consistency is achieved. 

The embedding method requires that we consider a rela- 
tively small number of surface layers; it is therefore a compu- 
tationally less expensive extension of DMFT in the presence 
of an interface compared to the slab method, in which the in- 
homogeneous system is simply represented as a finite number 
of layers. In this study, the number of surface layers is cho- 
sen to be N — 5 and we tested (by varying this number) that 
this number provides converged results. Our impurity solver 
is exact diagonalization, 16 where the bath is represented in 
terms of a finite number of levels, n s . For the case of phonon 
degrees of freedom we considered here, the infinite phonon 
space is also truncated allowing for a maximum number of 
excited phonons n p h- The typical values we considered for 
the bath level are n s = 8 and typical maximum number of 
phonons are n p h — 30 — 50. We tested that those numbers 
provide essentially converged results. For example changing 
n s from 8 to 9 changes z only by 4% for g = 0.5 which is 
close to the transition. For smaller g the error is smaller. 



III. RESULTS 

We use the technique explained in the previous section to 
study the Holstein model in a semi-infinite bipartite simple 
cubic lattice with in-plane translational symmetry and layer- 
dependent Hamiltonian parameters. We will work at half- 
filling (one electron per site), where any charge modulation 
is excluded by the particle-hole symmetry 17 and local occupa- 
tions on any layer, including the surface, coincide with the 



average filling, (n a ) = 1. We set the phonon frequency 
Slo = 0.2i, which puts the system in the adiabatic regime. In 
order to characterize the metal-insulator transition, we use the 
quasiparticle weight, z a = [1 — d£ Q (cj)/du;| w= o] _1 [S a (w) 
is the self-energy for layer a] whose vanishing marks the 
transition to the insulating state in which there is no spectral 
weight at the Fermi level. Another important quantity is the 
double occupancy, d a = (n a -fn a i), which is large in bipola- 
ronic states. 

We can model the inhomogeneity of the system and the 
different properties of the surface layer by introducing layer- 
dependent parameters. In particular, we can introduce differ- 
ent intralayer hopping t\\ or electron-phonon coupling g\ at 
the surface or we can tune the hopping between the surface 
and the second layer t\2. 

One immediately realizes that the actual behavior of the pa- 
rameters at the surface will depend on the specific properties 
of each material and on the geometry of the interface. On 
the other hand the aim of this work is to understand general 
tendencies of an electron-phonon system in the presence of a 
surface. Namely, we want to understand what happens when 
the surface is more metallic than the bulk and when the op- 
posite situation is realized. Therefore we will use one single 
parameter, t\\/t, to model the effect of all the others. The case 
tn < t will represent all the situations in which the surface 
is less metallic than the bulk, while tu > t will represent the 
opposite situation of a more metallic surface. 

In Fig. 2 we show the evolution of the layer-dependent 
quasiparticle weights, z a as a function of the ratio tujt for an 
e-ph coupling slightly smaller than the bulk critical coupling 
for the bipolaronic metal-insulator transition, g c ,bulk ~ 0.55. 

We first focus on the regime in which the surface is less 
metallic than the bulk, i.e., tu < t. We obviously find that 
the surface quasiparticle weight z a =i is smaller than that of 
the inner layers, although it does not vanish even for tu = 0, 
even though the e-ph coupling is very close to the bulk critical 
coupling. The finite (even if very small) value of z a —\ can 
be better identified by inspection of the self-energy which has 
to diverge to have a vanishing quasiparticle weight. It is seen 
from the inset of Fig. 2 that the surface self-energy S a= i (iui n ) 
goes up at low frequency and we expect extrapolation to zero 
as w n — > in a Fermi liquid manner while a divergence is 
excluded. 

This behavior is not unexpected because a metallic bulk is 
indeed able to determine an exponentially damped quasipar- 
ticle weight in the neighboring layers including the surface 
layer. This rules out the possibility to observe a truly in- 
sulating surface on top of a metallic bulk, even if polaronic 
effects will be amplified on the surface. 13 The double occu- 
pancy, shown in the bottom panel of Fig. 2, presents a strong 
enhancement at the surface layer with respect to all the other 
layers and the bulk, as expected by the reduced hopping which 
favors the e-ph coupling. The second and third layers present 
only small deviations with respect to the bulk. 

We now consider the case of a surface which is more metal- 
lic than the bulk, either because the surface e-ph coupling is 
smaller or as we now analyze, tu > t. The results, also re- 
ported in Fig. 2, show that, for very large tu z a —\ approaches 
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FIG. 2: Layer-dependent quasiparticle weight z a (top panel) and 
layer-dependent double occupancy d a (bottom panel) as a function 
of modified intra-layer surface hopping txi. We show results for the 
first three layers of the semi-infinite Holstein model with simple cu- 
bic (001) surface geometry and for the bulk quasiparticle weight. 
a = 1 indicates the topmost surface layer. The solid lines show 
bulk calculations. The insets show the imaginary part of the surface 
self-energy 7mEi(iw„) on the discrete mesh of the imaginary ener- 
gies ui n = (2n+l)n/p = 400) fortn/t = 0.1 andi n /t = 0.5. 
ImTi a= \{iuj n ) shows an upturn at small frequencies, compatible 
with a Fermi-liquid behavior. 



the free-electron value z a= \ = 1. This signals that the surface 
layer is essentially decoupled from the rest of the system and 
it supports an uncorrelated motion of the electrons. The rest 
of the system, however, remains strongly interacting and the 
a = 2 layer represents the new surface layer, the a = 3 layer 
becomes the first subsurface layer and so on. As is shown in 
Fig. 2, for all values of tn/t, the dependence of the quasipar- 
ticle weight in the subsurface layers on tn is comparatively 
weak and quickly diminishes with increasing distance from 
the surface. The behavior of the double occupancy confirms 
the decoupling of the topmost layer, which, for large values of 
tn/t, approaches the noninteracting value, 0.25. 

For g < g c .buik, the same qualitative behavior is observed 
by changing the inter-layer surface hopping ti2 ^ t or by 
changing the e-ph coupling at the surface gi ^ g. 

The decoupling between the surface and the bulk for large 
tn implies that one can in principle approach the situation 
where a metallic surface coexists with an insulating bulk. 
To investigate this possibility, we computed the quasiparticle 



FIG. 3: g dependence of quasiparticle weights z a of semi-infinite 
Holstein model for simple cubic lattice in the (001) orientation for 
enhanced intra-layer surface hopping. Surface transition at g — 
gc, surf ace- Bulk transition at g — g c 



weights as a function of g. In Fig. 3 we plot the quasiparti- 
cle weight vs. g for a moderately enhanced surface hopping 
rate, tn = 1.5t. Upon increasing g two different critical in- 
teractions are found. The first one marks the transition from 
a metallic to a bipolaronic insulating state at g c = g c ,buik, in 
which all the bulk quasiparticle weights (all layers except the 
surface) vanish. For larger e-ph interaction there is a range of 
values of g in which the bulk is a bipolaronic insulator while 
the surface is still metallic with a finite z a —i. Indeed, in this 
region some weight is induced in the subsurface layers. Since 
the low energy surface excitations cannot propagate into the 
bulk for g > g c ,buik and are instead reflected back to the sur- 
face for energies below the bulk excitation gap, the induced 
quasiparticle weight decreases exponentially with increasing 
distance from the surface. 

At a second critical coupling, g c , surf ace, the surface also 
becomes insulating and bipolaronic. For g > g c . sur f ace the 
entire system is in the bipolaronic insulating phase. A rather 
moderate enhancement of tn is sufficient to obtain a metallic 
surface phase. Obviously, a larger tn means that electrons 
in the first layer are more itinerant. A smaller surface coor- 
dination number clearly counteracts this mechanism. Conse- 
quently, we expect that a larger tn is needed to obtain a metal- 
lic surface state for more open surfaces, such as, for example, 
the (110) surface. The range of coupling g where a metallic 
surface coexists with an insulating bulk quickly increases as 
tn is increased. For tn — > oo the bulk energy scales become 
irrelevant and the electronic structure of the surface layer de- 
couples from the rest of the system. 

The overall results obtained here for a model with electron- 
phonon interactions are qualitatively similar to those of Ref . 8 
for a repulsive Hubbard model despite the fact that the nature 
of the transition is different in the two models. In both cases 
one can have an insulating surface coexisting with an insulat- 
ing bulk, which is a Mott insulator in the Hubbard model and 
a bipolaronic insulator in the Holstein model. Instead, in both 
models a more insulating surface gives rise to a single metal- 
insulator transition. Indeed the similarity between the two 
cases is not accidental. In the antiadiabatic limit the Holstein 
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FIG. 4: Phonon probability-distribution function for the first three 
layers of the semi-infinite Holstein model with nonuniform model 
parameters and the bulk phonon probability-distribution function for 
two different values of e-ph coupling strength, g. We used tn = 
1.5*. 



model becomes the attractive Hubbard model, which, at half- 
filling, can be mapped onto a repulsive model by a particle- 
hole transformation. Then the Mott transition of the repul- 
sive model becomes a pairing transition in which fermionic 
pairs are formed, 18 which in turn corresponds to the bipola- 
ronic transition in the Holstein model. Therefore the main dif- 
ference between the two models is the retarded nature of the 
electron-phonon interaction, as opposed to the instantaneous 
Hubbard interaction. In principle in the electron-phonon case, 
larger lattice distortions at the surface 13 could favor an insu- 
lating surface with respect to the case of the Hubbard model. 
However our results show that the dynamical nature of the in- 
teraction is not able to introduce qualitative differences with 
respect to a purely electronic model. 

To gain further insight about the region in which the surface 
remains metallic, we consider the behavior of the phonon dis- 
placement probability distribution function (PDF), P{x) = 
{(f>o\x) (x\(f>o) , where |c/>o) is the ground-state wave function 
and \x)(x\ is the projection operator on the subspace where 
the phonon displacement at a given site x has value x. 

This quantity is a measure of the distribution of the lo- 
cal distortions. 19 In the absence of e-ph interaction, P(x) 
is a Gaussian centered around x = 0. A small e-ph cou- 
pling slightly broadens the distribution which remains cen- 
tered around x = 0, implying that the coupling is not suffi- 



cient to give rise to a finite polarization of the lattice. Contin- 
uously increasing the interaction one eventually obtains a bi- 
modal distribution with two maxima at x = ±xrj. A bimodal 
PDF indicates that a certain number of the lattice sites are po- 
larized by the presence of electrons in such a way that the aver- 
age value of the polarization is greater than its fluctuations and 
therefore provide evidence for electron/phonon entanglement, 
i.e., a polaronic state. The point at which the phonon PDF be- 
comes bimodal is used as a marker of the polaron crossover, 20 
while P(x = 0) = can be used to characterize the transi- 
tion to the bipolaronic insulator (even if the vanishing of z is a 
more rigorous criterion). Figure 4 shows the phonon PDF for 
in = 1.5£ at two e-ph coupling values, one for g < g c ,buik 
(top panel) and one for g c ,bulk < 9 < 9c,surface (lower 
panel). In the first case the system is metallic but polarons 
are already formed in all the layers except the topmost one, 
which has a larger hopping amplitude, tu. Notice that the ap- 
pearance of polaronic distortion is not sufficient to make the 
bulk insulating, confirming that the polaron crossover and the 
bipolaronic transition do not coincide. Upon increasing the 
e-ph coupling and for g c ,bulk < g < g c ,surface, the phonon 
PDF of all layers except the surface go to zero at x = (bot- 
tom panel). This confirms the insulating phase of these lay- 
ers in this range of couplings while the surface layer shows a 
metallic state with polaronic character, as shown by the only 
slightly bimodal PDF. 



IV. CONCLUDING REMARKS 

We have investigated the effect of a surface on a strongly 
coupled electron-phonon system. We described this system 
with a Holstein model on a cubic lattice cut along one of the 
coordinate axes, assuming that the surface parameters are dif- 
ferent from the bulk ones. We used the ratio between the hop- 
ping within the surface layer and the hopping within bulk lay- 
ers, tn/t, to represent the effects of other nonuniform param- 
eters and considered both situations in which the surface is 
less metallic than the bulk (tn < t) and the opposite regime 
of a more metallic surface. Our focus is on the strong cou- 
pling regime, where we ask whether bulk and surface can be 
decoupled as far as the transport properties are concerned. In 
particular, for a given set of parameters, we ask whether one 
part of the system can be insulating while the other is metallic. 
We considered the system at half-filling, where the electron- 
phonon interaction can drive a bipolaronic phase transition. 
This choice inhibits charge transfer between the surface and 
the bulk. Our investigation, based on the embedding method 
for DMFT, 5 shows that a metallic surface can coexist with an 
insulating bulk when tn/t > 1 already for moderate values 
of this ratio while the opposite behavior is not realized even 
when the surface hopping is vanishing. The bulk excitations 
are always able to penetrate in the surface layer, even if they 
are strongly damped. Polaronic distortions, measured by the 
phonon distribution function, can be significantly different in 
the surface over a range of parameters. 

From the experimental point of view, even if the most typi- 
cal situation is that in which the surface is more insulating than 
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the bulk, evidence for a ferromagnetic metallic surface has 
been reported in antiferromagnetic insulating manganites 21 
and a surface insulator-to-metal transition has been observed 
in insulating NiS2- 22 In general these measurements are diffi- 
cult and rely on an interpretation of transport and magnetic 
measurements on samples with varying grain sizes. 23 Sim- 
ilarly, a use of photoemission and tunneling spectroscopies 
should be helpful to discern surface vs bulk properties. 
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